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Abstract: The relationships between a resistance curve (R-curve), the corresponding fracture 
process zone length, the shape of the traction/displacement softening law, and the propagation of 
fracture are examined in the context of the through-the-thickness fracture of composite laminates. 
A procedure that accounts for R-curve toughening mechanisms by superposing bilinear cohesive 
elements is proposed. Simple equations are developed for determining the separation of the 
critical energy’ release rates and the strengths that define the independent contributions of each 
bilinear softening law in the superposition. It is shown that the R-curve measured with a Compact 
Tension specimen test can be reproduced by superposing two bilinear softening laws. It is also 
shown that an accurate representation of the R-curve is essential for predicting the initiation and 
propagation of fracture in composite laminates. 
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1. Introduction 

To predict the propagation of damage in quasi-brittle materials such as composites, it is necessary 
to define damage evolution laws that account for the fracture energy dissipated in each damage 
mode. Thermodynamic consistency based on fracture toughness is necessary to ensure objectivity 
of the solution with respect to finite element mesh choice, to predict scale effects, and to determine 
the proper internal load redistributions. Most damage models, such as the Progressive Damage 
Model for Composites provided in Abaqus and typical cohesive elements, represent the evolution 
of damage with bilinear softening laws that are described by a maximum traction and a critical 
energy release rate. The shape of the softening law, e.g., bilinear or exponential, is assumed to be 
inconsequential for the prediction of fracture. 

The objective of the present work is to examine the relationships between the assumed shape of 
the traction/displacement softening law, the R-curve, its fracture process zone length, and their 
effects on the propagation of fracture. In Section 2, it is shown that material softening and R-curve 
are directly related to each other and that bilinear laws cannot accurately represent toughening 
mechanisms that cause an R-curve response. To address this difficulty, a procedure based on the 
superposition of bilinear softening laws is proposed that can account for all the mechanisms that 
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cause an increase in fracture toughness with crack growth. Finally, in Section 3 the fracture of a 
composite Compact Tension (CT) specimen characterized by a strong R-curve is analyzed with 
superposed cohesive elements. 


2. Softening Laws and the R-Curve 

In damage models such as the ABAQUS COH cohesive elements and the Abaqus Progressive 
Damage Model for Composites (Lapczyk, et ah, 2007), objectivity of the solution is achieved by 
representing a material’s softening response with bilinear traction displacement laws whose areas 
are equal to the critical energy release rate of the material in each mode of fracture, as shown in 
Figure 1. Each bilinear traction-displacement law is fully defined by three variables: an initial 
stiffness K p , a maximum strength <j c , and a critical energy release rate G c . 



The initiation of damage in a softening model occurs when the traction reaches the maximum 
nominal interfacial strength, o r . The integral of the tractions to complete separation yields the 
fracture energy release rate, G c . The length of the process zone (LPZ), l p , is defined as the distance 
from the crack tip to the point where the maximum cohesive traction is attained. Different models 
have been proposed in the literature to estimate the LPZ. Irwin, 1970, estimated the size of the 
plastic zone ahead of a crack in a ductile solid by considering the crack tip zone where the von 
Mises equivalent stress exceeds the tensile yield stress. Dugdale, 1960, estimated the size of the 
yield zone ahead of a mode 1 crack in a thin plate of an elastic-perfectly plastic solid by idealizing 
the plastic region as a narrow strip extending ahead of the crack tip that is loaded by the yield 
traction. Barenblatt, 1962, provided an analogue of the Dugdale plastic zone analysis for ideally 
brittle materials. Flillerborg, et al., 1976, introduced the concept of the LPZ as a characteristic of a 
material and estimated the LPZ for concrete to be 100 mm. Rice, 1979, estimated the LPZ as a 
function of the crack growth velocity. Bao, et al., 1992 and Cox, et al., 1994, examined the shape 
of the R-curve and LPZ for cracks bridged by different bridging laws. All of these models, 
regardless of whether they were developed for elastic or ductile conditions, predict that the LPZ 
has the form 


l P=y 
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where E is the Young’s Modulus in the transverse direction, and the nondimensional parameter y 
depends on the model and can have the values shown in Table 1. The LPZ represents an intrinsic 
characteristic of a material and the fact that it can range from a few nanometers for atomic bonds 
to several meters for concrete underlines the diversity in the response of different materials. 


Table 1. Value of parameter y in models for the length of process zone (LPZ). 
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In the fracture of concrete, where the LPZ can be as large as a few meters (Bazant, 2002) and 
where the scale from material characterization tests to a full-size dam usually spans three orders of 
magnitude, the inadequacy of linear elastic fracture mechanics (LEFM) to predict size effects was 
identified early on (ACI_Committee_224, 1992). It is now understood that for the analysis of con- 
crete structures one must take into account strain-softening due to distributed cracking, localiza- 
tion of cracking into larger fractures prior to failure, and bridging stresses at the fracture front. 

After examining the relationship between fracture testing of concrete and size effect, Bazant 
refined the concept of cohesive crack models where the contributions of small size and large size 
tractions is differentiated (Bazant, 2002). According to Bazant, a rather brittle initial mechanism is 
followed by a tougher secondary effect that acts at a lower stress level. As a result of extensive 
testing, the stress state of the initial fracture was determined to be about three times greater than 
the strength related to the secondary effect, while the toughness of concrete is approximately 2.5 
times the energy release rate associated with the initial fracture. 

Many common fracture processes in composites exhibit a fracture toughness that increases with 
crack growth - a response denoted as the resistance curve or the R-curve. In the presence of an R- 
curve, the toughness measured during crack propagation typically increases monotonically until 
the value stabilizes. In the case of delamination, the increase in toughness with crack growth is 
attributed mostly to fiber bridging across the delamination plane. Since it is generally assumed that 
fiber bridging only occurs in unidirectional test specimens and not in general laminates, the 
toughness of the material is taken as the initial toughness and the toughness for steady-state 
propagation is ignored (ASTM_D5528-01, 2002). 

In through-the-thickness fracture of composite laminates, the increase in fracture toughness with 
crack length is caused by a combination of damage mechanisms ahead of the crack tip and fiber 
bridging behind it. This R-curve makes it difficult to predict the effect of structural size on 
strength. Test results for laminated notched panels indicate that their strength cannot normally be 
predicted using a constant fracture toughness (MIL-HDBK-17-3, 2005). The degree to which scale 
effects in a composite laminate can be predicted by LEFM depends on its notch sensitivity, which 
is a property of both material and notch size, that can be expressed by the dimensionless ratio 77 of 
the notch dimension a over the LPZ: 

a [if 77 < 1 , notch - ductile 

77 = — { ( 2 ) 

/ [if 77 > 10 , notch - brittle 
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The load-carrying capability of notch-ductile components is dictated by strength; that of notch- 
brittle parts is dictated by toughness; and for the range in between these extremes, both strength 
and toughness play a role. The LPZ for polymeric composites is of the order of a few millimeters, 
so notched panels with notch lengths greater than approximately 1-10 cm are notch-brittle 
(Camanho, et al., 2007). Specimens that do not exhibit localized fracture planes, such as in the 
case of brooming, splitting, etc., may exhibit a much-enhanced notch-ductility. For instance, the 
notch-ductility is enhanced by multiple cracking and matrix splitting along the fiber direction. It 
can be shown that the fracture toughness of a composite laminate increases with fiber strength, Xf, 
and decreases with fiber/matrix shear strength, r, as (Bao, et al., 1992): 


G, oc - 


£/ 

Et 


( 3 ) 


The problem of length scales in bridged cracks has been examined by many authors, such as 
Foote, et al., 1986, Smith, 1989, Bao, et al., 1992, Cox, et al., 1994. In composite materials, more 
than one physical phenomenon is involved in the separation process, some acting at small opening 
displacements, which are confined to correspondingly small distances from the crack tip, and 
others acting at higher displacements, which will extend further into the crack wake. For example, 
when long-range friction effects are present, one might expect a law possessing a peak at low 
crack displacements corresponding to the tip process zone and a long tail at high crack 
displacements representing the friction that is transmitted through the layer of resin fragments in 
the wake of the crack (Yang, et al., 2005). 

Those processes whose length scale is negligible compared to the crack length can be represented 
accurately as a point of energy absorption at the crack tip in a LEFM formulation. According to 
Cox, et al., 1994, when the fracture toughness of a crack is due to a constant closing stress a b 
representing bridging plus a fracture toughness G T i p associated with the crack tip, the LPZ is 


/ 


p 


n EG b 
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where E is the transverse Young’s modulus and G b is the toughness resulting from the effect of 
the bridging tractions. When Gti p =0 , the expression for the LPZ reduces to a form similar Eq. 1: 


k EG b 


( 5 ) 


By comparing Equations 4 and 5, it can be seen that the presence of a toughness concentrated at 
the crack tip results in a process zone that is considerably shorter than one in which the bridging 
stresses resist the fracture alone. 


2.1 Relationships Between the Softening Law and the Shape of the R-curve 

The relationship between the shape of the strain-softening law and the R-curve has not received 
much attention. It is generally assumed that the fracture toughness of a composite is a material 
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constant that is independent of crack length. Toughening mechanisms such as fiber bridging that 
alter the characteristics of fracture propagation are often ignored. However, the assumption that G c 
is a constant is only a valid approximation when the LPZ is negligible compared to other 
dimensions such as the crack length. More strictly, the physics of stable crack growth should be 
viewed as the gradual development of a strain-softening bridging zone behind the crack tip that 
produces a stabilizing influence on cracking characterized by a raising crack growth resistance 
curve. In this way, “softening” may be a misnomer, for it is the softening mechanism that imparts 
“toughening” to the material’s fracture resistance. 

Foote, et al., 1986, developed an approximate expression for the crack-growth resistance curves of 
strain softening materials based on the assumption that the bridging stresses are insensitive to the 
crack opening profile and depend only on the distance from the crack tip. According to Foote, 
when the fracture toughness is the result of a toughness Gn P associated with the crack tip plus the 
toughening effect of linearly softening bridging stresses a h , the R-curve can be written as a 
function of the crack extension Acr. 


G R =G Tip+ l 


i^ajA 
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dx 


( 6 ) 


where S c is the critical displacement for the softening law and A a < l p . For the linear softening 
shown in Figure 1, 

A, =2^ 


( 7 ) 


By substituting Equations 1 and 7 into Equation 6 and performing the integration, we obtain an 
expression for the R-curve denoted as G p that is quadratic in Acr. 


G R {Aa) = 


G T i P +G b ^~ 


G Tip +G h 


\ 


2-^ 

V l P J 


for Act < I p 
for A a > l p 


( 8 ) 


where 1 p = y EG b / al . 

A simpler expression for the R-curve of a bridged crack was obtained by stress analysis by Bao, et 
al., 1992: 


c,M- 
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where l' p = y L EG b / a b . For linear softening, the LPZ coefficient y L is a function of crack 

extension, Aa. Bao calculated that y' ranges from 0.393 for short cracks to 0.732 for steady-state 
propagation. For a rectilinear (constant-stress) bridging law, the LPZ coefficient is a constant 
y i= 0.393, as illustrated in Figure 2. Given the same <7* and G b , the LPZs for linear and uniform 
bridging stresses in Bao’s model differ by approximately a factor of two. 



Figure 2. R-curves for rectilinear bridging and softening (after Bao, et al., 1992). 

By equating Equations 8 and 9, the following relation between yand y L is obtained: 


y L (Aa) 


Y 



(10) 


where y is a constant and y L is a function of Aa. In the present work, we linearize Equation 9 by 
assuming that y L is a constant, regardless of the softening law. Assuming that the mean value of 
y L is reached at Aa=l /2, Equation 10 gives y L = fy and l p = -|/ . The corresponding 
linearized R-curve for linear softening is shown in Figure 2. 


2.2 R-curve by Superposition of Two Bilinear Softening Laws 

The preceding section indicates that when l p «a, all bridging mechanisms can be assumed to be 
acting at the crack tip and the R-curve effect can be ignored. Therefore, for a vast number of 
problems, the shape of the softening law has a negligible influence on the prediction of crack 
propagation. Flowever, in problems where the notch sensitivity /] defined by Equation 2 is lower, 
the increase in the fracture toughness with crack growth R-curve effect becomes important and 
more complicated softening laws may be necessary to capture the correct bridging response. 

When softening laws more complex than bilinear are necessary, new damage models have to be 
formulated. Alternatively, multilinear softening laws can be obtained by combining two or more 
bilinear softening laws, as illustrated in Figure 3. The two underlying bilinear responses may be 


6 


2008 ABAQUS Users ’ Conference 


seen as representing different phenomena, such as a quasi-brittle delamination fracture 
characterized by a short critical opening displacement <5 cl , combined with fiber bridging, 
characterized by a lower peak stress, a longer critical opening displacement <? c 2 , and a 
correspondingly shallower R-curve. 

To obtain a desired experimentally-determined R-curve from the combination of multiple bilinear 
laws, it is necessary to determine the proportions of the contributions of each bilinear law to the 
total critical energy release rate and the strength. For the superposition shown in Figure 3, we 
define the proportions a cX =no c , a cl =(1 -n)cr c , G ] =mG c and G 2 = (1 -m)G c with 
0 < n, 7« < 1 , so that 

G c =G l +G 2 and a c = cr cl + <j c2 (11) 



Figure. 3. Trilinear softening law obtained by superposing two bilinear laws. 

The objective of the present work is to develop a procedure for determining the stress ratio n and 
the toughness ratio m that approximate an experimentally-determined R-curve. Using Equation 8 
as the basis for the R-curve, we define the following empirical form of the R-curve that results 
from the sum of two bilinear softening laws: 


?r(A a) = nG c — 


A a 
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2 n An 

777 L 


+ (1-77 )G C ^~ 
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(l — 77 ) Afl 
(1 - 777) l c 


=G\ if A a> — l c 


=Gt if Aa >— l 
1 -n 


( 12 ) 


where l c =y EG C / cr? is the LPZ for bilinear softening. It can be seen that Equation 12 satisfies 
the following conditions: 

• if 77=777 (which is equivalent to having - &C2 in Figure 3), then Equation 12 is 
independent of n and m and it reverts to Equation 8. 

• when A a =l c min, the first term on the RHS of Equation 12 is equal to Gj. 
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when A a =/ c ( 1 - »;)/(! - n ) , the second term on the RHS of Equation 12 is equal to G 2 . 


If the softening laws 1 and 2 are ordered such that m / n < (1 - m) /(I - n ) , then the LPZ for trilinear 
softening, defined as the crack extension necessary to reach the steady-state value of Gr=G c is 


jNL _ 1 — . Otheiwise, l NL = —l 
p l-n p n 


(13) 


Consequently, /* > l c , i.e., the length of the process zone for trilinear softening cannot be shorter 
than that for bilinear softening with the same critical energy release rate (ERR) G c . 

A simpler expression for the R-curve for the sum of two bilinear softening laws can be constructed 
using Equation 9 and the approximation y ~\y. The proposed R-curve has the form 


Gr(Aci) = MIN 
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As with Equation 12, three separate ranges can be identified. If the softening laws 1 and 2 are 
ordered such that m/n<( 1 - m) /(I - n ) , then: 


c«(a«)= 


3 A a 

Cr r 

2 L c 


3 A a 


2 L 


if A a<^ — l c 

n 


Gi+(l-n)-—G c if 


2 m 2 / 

4 — L < A a< 4 L 


n 


3 


l-n 


if A 

3 l-n c 


(15) 


Therefore, the linearized model of the LPZ for trilinear softening is 

,L _2 l ~ m , 

p 3 l-n c 


(16) 


In the remainder of the present work, the LPZ coefficient of Rice’s model, y=0.884, was chosen 
based on the results of the cohesive element analyses of Turon, et ah, 2007. 


3. Compact Tension Specimen: Test and Analysis 

To obtain the toughness of fiber fracture in T300/913 graphite/epoxy, Pinho, et ah, 2006, used the 
Compact Tension specimen shown in Figure 4 with a stacking sequence of [0/90] 8s . By examining 
the micrographs of the fracture zone, Pinho observed that the mode of fracture of the fibers 
perpendicular to the fracture plane evolves from pure fiber fracture to a combination of fiber 
fracture, fiber pullout and fiber bridging during crack propagation. By performing a J-integral 
calculation of the energy release rate (ERR) as a function of the observed crack tip position, Pinho 
determined that the ERR for fracture initiation is about half of the ERR for steady-state 
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propagation and that it takes about 1 1 mm. of crack growth to reach the steady state propagation. 
Pinho’s test results show that the crack propagates in short bursts, indicating that local dynamic 
effects contribute to the crack propagation. Therefore, it is the upper bound of the ERR curve that 
most closely reflects the quasi-static fracture toughness of the laminate. This upper bound is 
approximately G = 1 80 kJ/m 2 . 

3.1 Analysis of Fracture of CT Specimen with Cohesive Elements 


To investigate the effects of different softening laws on the propagation of crack in Pinho’s CT 
tests, the shell model shown in Figure 5 was created. Symmetry is assumed along the plane of 
crack propagation. A row of four-node cohesive elements (Turon, et ah, 2006) is used to model 
softening and crack propagation. The length of the elements along the fracture line is 0.208 mm. 
The material properties listed by Pinho were used: £,, = 131.7 GPa. £\,=8.8GPa, G /2 =4.6GPa, 
and u / 2 =0.32. The corresponding Young’s Modulus for the [0/90] gs laminate is £=70.62 GPa. The 
value of the longitudinal strength is taken from Camanho, et al., 2007 as cr,=2,326 MPa. 



Figure 4. Compact Tension Figure 5. Shell model of Compact Tension 

Specimen (after Pinho, et al., 2006). Specimen with cohesive elements for 

fracture. 


To improve the convergence rate of the load incrementation procedure, the analysis was run 
dynamically using Abaqus/Std. Two levels of Rayleigh stiffness-proportional damping were 
applied: a low value of j3 R = 2.5 10' 4 and a moderate value of J3 R = 2.5 10" 3 . The initial stiffness 
predicted by the model was found to be 5% higher than the results reported by Pinho and this 
difference was attributed to the compliance of the test machine. Therefore, the compliance was 
estimated and subtracted from the published displacements as 

^corr — ^ meas ^'mach^A ( ^ ^ 

where F A is the applied load and the machine compliance is estimated to be C mach =0.13 mm / kN ■ 

Because the model is symmetric, only half of the experimental critical steady-state ERR G, = 180 
kJ/m 2 is applied in the definition of the cohesive elements. The predicted results are shown in 
Figure 6. It can be observed that the results using the higher damping coefficient correspond to an 
upper bound of the results with the reduced damping, which indicates that dynamic effects play a 
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role in the jagged force plot. Also, it can be observed that the predicted strength for both damping 
levels is 29% higher than the experiment. However, for applied displacements larger than 2.4 mm, 
the analyses correlate well with the experimental results. 


The crack tip position, a, extracted from the analysis with low damping was recorded as a function 
of the model stiffness, K A =F A /S A . The results are plotted as {k a ) x vs. a, where the exponent y 
is chosen to make the curve a straight line, as shown in Figure 7. The curve fit of the crack length 
as a function of the compliance is then 


a = 


_(K a Y-£_ 


a 


(18) 


where the best fit coefficients were determined to be a=-2.717, /?=137.87, and j=0.53. 



Applied Displacement, mm. 
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- y = 1 37.87 - 2.71 69x R= 0.99994 
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35 40 45 

Crack Length, a, mm. 


Figure 6. Force-displacement Figure 7. Plot of FEM model stiffness 

response of CT specimen. K A as a function of crack length a. 


The R-curve for the analysis can be obtained by either the Area Method (AM) or the Modified 
Compliance Calibration method (MCC) (ASTM D5528-01, 2002). With the AM method, the 
instantaneous toughness of the material is obtained by dividing the energy released between two 
consecutive points in the force-displacement curve by the new fracture surface created. The 
expression is: 


G f+1 

R 


( d i+i F i- d i F i+\) 

2 -t 0 •(«,■+!-«/) 


(19) 


where t" is the thickness of the 0-degree plies and the crack length increment (a i+ i~ a,) is 
calculated using Equation 18. The toughness of the 90-degree plies is three orders of magnitude 
lower than that of the fibers, so it is neglected in the present calculation. 

The equation for the energy release rate with the MCC method is 
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where the rate of change of the compliance C with crack length a can be obtained by 
differentiation of Equation 18, which gives for the MCC method: 


G K (a) = 


F% a(j3 + aa) ( 1+ /z) 
2 • t° I 


(21) 


where the crack length a can be obtained from Equation 18. The R-curves calculated from the 
finite element analysis results using the AM and MCC methods for a damping ratio of f) R = 2.5 10" 3 
are shown in Figure 8. These results indicate that G R rises rapidly from zero and quickly settles 
into a plateau of steady-state propagation where Gr = 180 kJ/m 2 . A comparison of the R-curves 
obtained with the MCC and the AM methods indicates that the two methods provide nearly 
identical results. Elowever, the MCC method is evaluated point-wise rather than through backward 
differentiation, so it is potentially more accurate than the AM method. Since the MCC method is 
also easier to apply than the AM method, the MCC method is used in the remainder of this work. 

The R-curves predicted using the nonlinear model G R L given by Equation 8 and the linear model 
G l r given by Equation 9 with y=0.884 and G T i P = 0 are also shown in Figure 8. The LPZ is l p = 2.08 
mm. It can be observed that the R-curves predicted with the two approximate formulae correlate 
well with each other. 

The R-curves obtained with the MCC method for damping ratios of /3 R = 2.5 10’ 4 and j3 R = 2.5 10" 3 
are shown in Figure 9. It can be observed that there is a significant amount of scatter in the lightly 
damped results and that the results of the analysis with higher damping provide an approximate 
upper bound to the results with less damping. The non-monotonic appearance of the crack 
propagation is due to the dynamic forces that perturb the instantaneous stiffness that is used in the 
calculation of the crack length (Equation 18). 

The MCC method was also applied to the force displacement test results of Pinho with the 
compliance calibration coefficients determined by the present finite element analysis. The results 
are shown in Figure 10. The R-curve initially rises rapidly from zero to about 100 kN/m 2 and then 
continues to rise at a slower rate until steady-state propagation at 180 kN/m 2 . The LPZ for steady- 
state propagation is about 1 1 mm. After 20 mm of propagation, the data should be truncated 
because the proximity of the edge induces an artificial reduction in the measured toughness G R . 

A comparison of the predicted R-curves in Figure 8 and the experimental results in Figure 10 
indicates that the finite element analysis with a single bilinear cohesive element cannot reproduce 
the characteristics of the measured R-curve. The fact that the rising portion of the R-curve is 
severely overestimated by the analyses explains why the analyses also overpredict the initiation of 
fracture as shown in Figure 6. However, the models do predict the correct critical ERR for steady 
state propagation which allows the steady-state load-deflection response to be predicted as well. 

In the next section, superposed cohesive elements are used to obtain a better representation of the 
R-curve for a more accurate prediction of the strength of the CT specimen. 
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Crack Length, a, mm. 


Figure 8. R-curve G R obtained from 
FEM using MCC and AM methods. 


Figure 9. Effect of Raleigh damping (3 R 
ratio on R-curve G R . 



Figure 10. R-curve obtained by MCC using Pinho’s test data. 

3.2 Fracture Analysis with Superposed Cohesive Elements 

The values of n and m define the properties of the two superposed cohesive elements that must be 
used to reproduce an experimental R-curve can be obtained graphically using the linear R-curve 
given by Equations 15 and 16. The outline of the experimental R-curve in Figure 10 can be 
idealized as being composed of a critical ERR associated with the processes occurring in the 
vicinity of the crack tip, G/=100 kJ/m 2 , plus a critical ERR associated with a bridging mechanism, 
G 2 = 80 kJ/m 2 . Therefore, at least two cohesive elements are necessary to represent the cohesive 
tractions where m=Gi/G c =0.556. The strength ratio n is determined by solving Equation 16: 
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For the material properties of T300/913 presented in Section 3.1 and an LPZ =11 mm, the 
stress ratio is n= 0.944. If Hillerborg’s LPZ factor /=1.0 is used instead of 0.884, the stress ratio 
becomes n =0.937. Therefore, the LPZ factor has a small influence on the predicted stress ratio. 

The values for n and m can also be obtained with a more formal curve fitting of the experimental 
R-curve using the nonlinear model defined by Equations 12 and 13. Flowever, the linear model is 
simpler to apply and it fits the experimental data well. Therefore, the nonlinear model was not 
used in the present analysis. 

To verify that the model with superposed cohesive elements can reproduce the correct R-curve, 
the finite element model of the previous section was modified with doubled cohesive elements. No 
additional modeling effort is required since the superposition is done with a single command line: 

*ELCOPY , ELEMENT SHIFT=<SHIFT_N>, OLD SET=COHl , SHIFT NODES=0 , NEW SET=COH2 

which creates a new set of cohesive elements without creating new nodes. 

The R-curve extracted with the MCC method from the finite element analysis is shown in Figure 
11. The R-curves Gr L and G' r predicted by Equations 12 and 14 are also shown in Figure 11. 
Using Equations 13 and 16, the LPZ for the nonlinear model is lp L =-§-11 rnm = 16.5 mm. It can be 
observed that in spite of the difference in the LPZ between the two models, the curves described 
by the two equations correlate well with the R-curve extracted from the finite element analysis and 
with the outer envelope of the test results. 

After more than 15 mm of propagation, some divergence occurs between the FEM and the test 
results, which is due the proximity of the crack tip from the edge of the specimen. The specimen is 
fully fractured when a is equal to 51 mm (not shown in the figure). 




Applied Displacement, mm. 


Figure 11. R-curves for analytical 
models, FEM, and experiment. 


Figure 12. Load-displacement response 
of CT specimen calculated using 
superposed cohesive elements. 
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The load-displacement response curves of the CT specimen predicted with two bilinear cohesive 
elements with n = 0.944 and m = 0.556 are shown in Figure 12 for damping ratios J3 R = 2.0 10" 4 and 
jBr = 2.5 10"'. By comparing these results to those shown in Figure 6 for a bilinear cohesive law, it 
can be observed that the use of superposed elements has reduced the error in the predicted strength 
of the CT specimen from 29% to 2.8%. This dramatic improvement in accuracy was obtained 
without modifying any material property used in the analysis and without any additional modeling 
effort beyond that required to find the strength ratio n and the toughness ratio m for the separation 
of the fracture response into two parts. 

4. Concluding Remarks 

The importance of the R-curve, and in particular the length of the process zone, on the prediction 
of through-the-thickness fracture of composite laminates was examined. The fracture of a 
Compact Tension specimen used by Pinho to measure the fracture toughness of carbon fibers was 
analyzed using cohesive elements. It was found that a single bilinear softening law is insufficient 
to account for the multiple damage mechanisms and the toughening due to fiber bridging and fiber 
pullout that occurs during the fracture of a composite laminate. However, it was shown that by 
combining two bilinear cohesive elements of different properties, it is possible to represent the 
resistance curve that is necessary for predicting the initiation and propagation of fracture. 

Two formulae, one linear and the other nonlinear, that predict the shape of the R-curve that results 
from the superposition of two bilinear softening laws are proposed. By fitting these formulae to a 
measured R-curve, it is simple to determine the parameters that define the two bilinear softening 
laws. This fitting can be performed by visual inspection of the R-curve for the linear law, or it can 
be done by a more formal optimization in the case of the nonlinear law. It was shown that with the 
use of superposed cohesive elements the error in the predicted strength of a CT specimen is 
reduced from 29% to 2.8% without modifying any material property and without any additional 
modeling effort. In addition, the proposed approach relies on available bilinear cohesive elements 
so additional finite element development is unnecessary. 
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